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The three dimensional ferromagnetic Ising model with a random field (RFIM) shows a 
phase transition to long range order at a critical temperature for small enough field strength 
||T||. However, the nature of this transition is still unclear; even the question of whether it 
is first or second '^j^ order remains unsettled. The droplet theory of Villain ^ and 
Fisher |0 (see also Bray and Moore J^) develops a self-consistent picture of the transition 
as well as a set of scaling relations among the critical exponents. Existing numerical studies 
have been unable to test the validity of these scaling relations because not all the exponents 
were calculated for any of the relations. The aim of this paper is to determine all critical 
exponents within a single numerical simulation in order to test the scaling relations predicted 
by the droplet picture. 

The droplet picture makes other predictions which are relevant to our simulations. One 
is that at and below the transition temperature T^, the susceptibility is expected to have 
large sample to sample fiuctuations 0. We therefore need to average over a large number 
of samples to get a good statistics. Another prediction is thermally activated dynamical 
scaling |^,^ resulting in a dramatic slowing down in the critical region. This means very 
long equilibration times. For these reasons we had to confine ourself to modest lattice sizes 
and the critical exponents will be obtained via finite size scaling. 

The Hamiltonian of the system is given by 

n = -J2s^s,-Y.h,s,, (1) 

(ij) i 

where Si = ±1 are Ising spins and the first sum runs over all nearest neighbor pairs on an 
L X L X L simple cubic lattice with periodic boundary conditions. The random fields hi 
in the second sum, running over all sites, take random values with the discrete probability 
distribution 

Pih,) = ^5ihi - K) + U{hi + K) . (2) 

The Monte Carlo (MC) simulations were performed on a transputer array using 40 T414 
transputers. We were able to obtain high performance by using a multi-spin coding algo- 
rithm described in ||T^ in which each transputer simulates 32 physically different systems 



in parallel, each with a different random field realization. This is somewhat different from 
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the implementation of multi-spin coding which was applied to the RFIM in [jTT|. For each 
run at fixed temperature, T, and field-strength, hj., we performed a disorder average over 
1280 samples. An average over such a large number of samples is necessary because the 
susceptibility is highly non self-averaging, as mentioned before. The simulations were done 
for fixed ratio hr/T at different temperatures. 

To check equilibration we simulated two replicas of the system: one starting from an 
initial configuration with all spins up and one with all spins down. We assumed that the 
system has reached equilibrium when the magnetization measured for both replicas is the 
same (within the errorbars). The time needed for equilibration of all 1280 systems varied 
much with the system size and temperature — in case of the largest size {L = 16) we used 
up to 0.5-10^ MC-steps for equilibration and 1.5-10^ MC-steps for measurements. All of the 
samples were equilibrated for L < 10. For larger sizes the number of nonequilibrated samples 
generally varied between 1% and 3%. The contribution of these samples was estimated to 
be less than the error bars in the points so so no significant error was made by including 
them. The only exception to this was for L = 16, hr/T = 0.5 for which 5% of the samples 
were not equilibrated which gave a significant error in the susceptibility, though not for the 
other quantities. We therefore ignored this data point when analyzing the susceptibilty. 

For each sample and each replica {a,b) we recorded the average magnetization per spin 
(Ma^h), its square (M^^,), the average energy per spin (-Ea.b) and its square {E^fj). The 
angular brackets, (...), denote a thermal average for a single random field configuration. 
From these data we get the specific heat per spin, C, the susceptibihty Xi the disconnected 
susceptibility, Xdis, and the order parameter, m, as follows: 

[C],. = n{[{E%^^ - m\^]/T^ , [m],. = [|(M)|],. , 

(3) 

[X]av = iV{[(M2)],, - [(M)2],,}/T , [Xd.]av = N[{M)\^ , 

where [. . .]av denotes the average over different random field configurations. 

The procedure we used to extract the critical exponents is the following: Let t = T — T^ 
(the deviation from the critical temperature), then the finite size scaling functions for the 
above quantities read: 
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(4) 

T [xU = L'-^x{tL'^n , [Xdis]av = L'^-'^XdUtL'/n , 
where a is the specific heat exponent, /3 the order parameter exponent, u the correlation 
length exponent and t] and rj describe the power law decay of the connected and disconnected 
correlation functions, see e.g. [|l|. Note that the susceptibility exponent is given by 7 = 
(2 — ?7)z/. The scaling function C{x) has a maximum at some value x = x*. For each lattice- 
size we estimate the temperature T*(L), where T^[C]av is maximal. Since t*(L) L^^" = x* we 
obtain in this way the critical temperature Tc and the correlation length exponent u from: 

t*{L) = T*{L) -T, = x*L-^/^ . (5) 

We denote the value of T*{LYC at this temperature T*{L) by C* and similarly for 
the other quantities in Eq. (^. In the vicinity of x* the scaling function C{x) can be 
approximated by a parabola. Therefore three temperatures near the maximum of the specific 
heat are enough to determine the values of T*(L) as well as [C*]av etc. Our results for the 
exponents obtained in this way are summarized in Table 1. For illustration, we show the 
results for T*{L), [C%,, [x*U, [m%, and [xli.U for hr/T = 0.35 in Fig. la-d. Several 
comments have to be made: 

1) The higher the field strength the harder it is to equilibrate the samples. However, 
the lower the field the less pronounced is the random field behavior for small lattice sizes 
because of crossover from pure Ising model behavior. Therefore the investigation of larger as 
well as smaller ratios h^/T did not seem to be advisable to us. If the transition is of second 
order and no tricritical point occurs along the critical line (Tc, he) the exponents should be 
universal, i.e. independent of the value of hr/T. 

2) The shift of T*{L) with respect to Tc becomes smaller for low field strength, so it is 
harder to determine the exponent u. In case of h^/T = 0.25 it was not possible to perform 
an acceptable fit for T*{L) according to equation (^. The values of u obtained for the other 
ratios hr/T are somewhat higher than that obtained in |Q, where z/ = 1.0 ± 0.1. 

3) We did not find any indication of a divergence, even logarithmic, of the specific 
heat, so a is negative. This is different from what is found experimentally |1|,§|, where 
the specific heat diverges logarithmically, corresponding to a = 0. Furthermore, in our 
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simulations a seems to get more negative with increasing ratio h,./T . This may indicate 
that it is difficult to determine a when a is negative because non-singular (but temperature 
dependent) background terms can give a significant contribution to the specific heat. 

4) The order parameter [m*]av shows only a very small size dependence, and does not 
approach zero but limL_,oo["^*]av ~ 0.52, 0.50 and 0.47 for hr/T = 0.5, 0.35 and 0.25, 
respectively (see the inset of fig. Id). This indicates that (3 = and that the transition is 
first order. This seems to contradict our results for the specific heat since the specific heat is 
expected to diverge as L"^ JTSf at a first order transition, because of the latent heat, whereas 
our specific heat data seem to saturate for large L. Perhaps the coefficient of is zero 
(though we see no symmetry reason for this) or is so small that L'^ behavior would only be 
seen for larger sizes. 

5) For the exponent rj we get a best estimate that is slightly higher than 1/2, which is the 
value obtained below and also the value at if the transition is first order . However, 
the value t] = 0.5 is not excluded by our data. For hr/T = 0.5 we had to exclude the size 
L = 16 from the analysis since 5% of our samples were not equilibrated and the contribution 
of these samples was larger than the error bars. Our estimates for t] are consistent with that 
obtained in rj = 0.5 ± 0.1. 

6) The exponent rJ for the disconnected susceptibility turns out to be equal to one, so 
that the scaling relation 

P= ^d-4 + rj)u/2 (6) 
is fulfilled (as indicated in the table). We also see that the Schwartz-Soffer [T^ inequality 



V<2rj, (7) 

holds as an equality within the error bars. In it was found that rJ = 1.1 ± 0.1. 

7) In the droplet picture [^,0 6 = 2— rj+r^ is called the violation of hyperscaling exponent. 
The hyperscaling relations then have the spatial dimension, d, replaced hj d — 6, e.g. 

2-a = {d-9)u. (8) 

As indicated in the table this equality seems not to be fulfilled though the error bars are 
quite large and our estimate for a might be affected by temperature dependent background 



terms, as discussed above. The estimates of both sides of this equation cannot be made 
without knowing both a and u but for h^/T = 0.25 we only determined the ratio. The 
entries in the table for 2 — a and {d — 6')z/ are therefore left blank for hr/T = 0.25. 

8) One of the main predictions of the droplet picture is a long tail in the distribution 
of the susceptibility x fo^' samples of size L aX T = 0. An analysis of this distribution 
extracted from our results for the 1280 samples confirms the existence of this long tail. 
Figure ^ shows the histograms for the probability distribution P{x) close to the temperature 
T* (T=3.80 for L=8 and T=3.75 for L=16) for /;,j,/T=0.35. The second moment of this 
distribution [x*^]av; shown in the inset of fig. Ic, scales like L'', with = 3.8 ±0.1 (for 
hr/T = 0.35), which is larger than the square of the mean L^"^'' ~ L^-^, but somewhat 
smaller than the predicted value ( = 6 — fj — t] ^ 4.4 0. We attribute this difference in ( to 
the number of samples being too small to catch a sufficient number of rare samples which 
dominate the higher moments. 

To conclude, while the data for the magnetization and disconnected susceptibility in- 
dicate a first order transition fairly convincingly, the specific heat seems to saturate to a 
finite value so there is no detectable latent heat. It is interesting to ask if the order of the 
transition might be different for a different random field distribution, since mean field theory 
predicts [0 that the transition becomes first order for large fields for the ±h distribution, 
but not, for example, for the Gaussian distribution. Since the multi-spin coding technique 
that we used does not work for a continuous distribution of fields, the answer to this question 
needs an even larger computing effort. Nevertheless we are currently attempting to carry 
out similar calculations for the Gaussian distribution. Our results are consistent with the 
Schwartz-Soffer inequality, Eq. (|^), being satisfied as an equality, and support the scaling 
relation, Eq. @. The scaling relation involving the specific heat, Eq. (P), does not seem to 
be satisfied, though our values for a may only be effective exponents particularly since we 
find a is negative and so a more detailed determination of non-singular background terms 
might be necessary to determine a accurately. Our results do support the prediction of the 
droplet theory that there are large sample to sample variations in the susceptibility at T^. 
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FIGURES 



FIG. 1. The results of least square fits to data obtained by the procedure described in the text 
for hr/T = 0.35. The points indicated by diamonds (o) correspond to lattice sizes L=4, 6, 8, 10, 
12, 16 — from right to left in (a) and (b) and left to right in (c) and (d). With the exception 
of the susceptibility we also inserted data for L = 24 with squares (□), which we obtained by 
using a the same algorithm but on a CRAY Y-MP instead of the transputer array and which are 
averaged over only 64 samples. The L = 24 data were not used for the least square fits, (a) The 
temperature T*(L) of the specific heat maximum versus L~^l^ with v = 1.64 and = 3.552, see 
Eq. (jsp. (b) Specific heat [C*]av = L'^/''C{x*) versus L'^/" with a = 1.04 and v as in (a), see Eqs. 
@ and (js]). The specific heat appears to saturate for L ^ oo at a value of limL_,oo[C'*]av — 25.3. 
(c) Susceptibility [x*]av = L'^~^x{x*) in a log-log plot. The slope of the straight line is 2 — ry 
with rj = 0.53. The inset shows the second moment [x*^]av of the probability distribution P{x) 
in a log-log plot. The slope of the straight line is C = 3.82. (d) Disconnected susceptibility 
[Xdislav = L'^'^'^Xdisix*) in a log-log plot. The slope is 4 — rJ with rJ = 1.0. The insert shows the 
magnetization [m*]av = L~^/'^fh{x*) as a function of a L (the scale of the x-axis is logarithmic). 
The straight line is the extrapolation to [m*]av(-Z^ = oo), which is clearly nonzero and so /? = 0. 

FIG. 2. The histograms for the probability distribution P{x) of the susceptibility for different 
L = 8 and 16 with hr/T = 0.35. The temperatures are chosen to be as close to T*{L) as possible: 
T = 3.80 for L=8 and T = 3.75 for L=16. 

The y-axes of the inserts are scaled differently to emphasize the long tail of the distribution. This 
feature originates in the rare samples with the extremely large values of the susceptibility scaling 
with the volume of the system (since 4 — 77 3) . 
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TABLES 

TABLE I. The critical exponents obtained via finite size scaling according to the procedure 



described in the text. 



hr/T 


0.25 


0.35 


0.5 


Tc 

V 

a. 

V 
V 

Q = 2-r} + ri 


3.9 ±0.1 

|^=-0.50 ±0.05 

0.60 ±0.03 
0.97 ±0.08 
1.6 ±0.1 


3.55 ±0.05 
1.6 ±0.3 
-1.0 ±0.3 
0.56 ±0.03 
1.00 ±0.06 
1.6 ±0.1 


3.05 ±0.05 
1.4 ±0.2 

-1.5 ±0.3 
0.6 ±0.1 
1.04 ±0.08 

1.6 ±0.1 


(d-4 + 7?)z^/2 




±0.05 




±0.05 




±0.05 


{d~e)u 

2-a 




2.3 ±0.5 
3.0 ±0.3 


2.0 ±0.6 
3.5 ±0.3 
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Critical Temperature 
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Fig. la 
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Fig. lb 
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Fig. Ic 
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Fig. Id 
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